﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using SpiceNet.Circuits;
using SpiceNet.Diagnostics;
using System.Numerics;
using SpiceNet.Simulations;
using SpiceNet.Integration;
using SpiceNet.Parameters;

namespace SpiceNet.Models
{
    /// <summary>
    /// An instance of a diode
    /// </summary>
    public class Diode : Instance
    {
        /// <summary>
        /// Get parameters
        /// </summary>
        public override Dictionary<string, IP> ParamTable => _pTable;

        /// <summary>
        /// Private variables
        /// </summary>
        private int DIOposNode, DIOnegNode, DIOposPrimeNode, DIOstate;
        private double DIOcap;   /* stores the diode capacitance */
        private double DIOtJctPot;  /* temperature adjusted junction potential */
        private double DIOtJctCap;  /* temperature adjusted junction capacitance */
        private double DIOtDepCap;  /* temperature adjusted transition point in */
        /* the curve matching (Fc * Vj ) */
        private double DIOtSatCur;  /* temperature adjusted saturation current */
        private double DIOtVcrit;   /* temperature adjusted V crit */
        private double DIOtF1;      /* temperature adjusted f1 */
        private double DIOtBrkdwnV; /* temperature adjusted breakdown voltage */
        private Noise.InstanceNoiseData[] DIOnVar = new Noise.InstanceNoiseData[4];

        #region Parameters
        private enum _c
        {
            DIO_AREA = 1,
            DIO_IC = 2,
            DIO_OFF = 3,
            DIO_CURRENT = 4,
            DIO_VOLTAGE = 5,
            DIO_CHARGE = 6,
            DIO_CAPCUR = 7,
            DIO_CONDUCT = 8,
            DIO_AREA_SENS = 9,
            DIO_POWER = 10,
            DIO_TEMP = 11,
            DIO_QUEST_SENS_REAL = 12,
            DIO_QUEST_SENS_IMAG = 13,
            DIO_QUEST_SENS_MAG = 14,
            DIO_QUEST_SENS_PH = 15,
            DIO_QUEST_SENS_CPLX = 16,
            DIO_QUEST_SENS_DC = 17,
            DIO_CAP = 18
        }
        private static Dictionary<string, IP> _pTable = new Dictionary<string, IP>
        {
            { "off", new IP(IF.IOPU, (int)_c.DIO_OFF, T.FLAG, "Initially off") },
            { "temp", new IP(IF.IOPU, (int)_c.DIO_TEMP, T.REAL, "Instance temperature") },
            { "ic", new IP(IF.IOPAU, (int)_c.DIO_IC, T.REAL, "Initial device voltage") },
            { "area", new IP(IF.IOPU, (int)_c.DIO_AREA, T.REAL, "Area factor") },
            { "sens_area", new IP(IF.IP, (int)_c.DIO_AREA_SENS, T.FLAG, "flag to request sensitivity WRT area") },
            { "vd", new IP(IF.OP, (int)_c.DIO_VOLTAGE, T.REAL, "Diode voltage") },
            { "id", new IP(IF.OP, (int)_c.DIO_CURRENT, T.REAL, "Diode current") },
            { "c", new IP(IF.OPR, (int)_c.DIO_CURRENT, T.REAL, "Diode current") },
            { "gd", new IP(IF.OP, (int)_c.DIO_CONDUCT, T.REAL, "Diode conductance") },
            { "cd", new IP(IF.OP, (int)_c.DIO_CAP, T.REAL, "Diode capacitance") },
            { "charge", new IP(IF.OPU, (int)_c.DIO_CHARGE, T.REAL, "Diode capacitor charge") },
            { "capcur", new IP(IF.OPU, (int)_c.DIO_CAPCUR, T.REAL, "Diode capacitor current") },
            { "p", new IP(IF.OPU, (int)_c.DIO_POWER, T.REAL, "Diode power") },
            { "sens_dc", new IP(IF.OPU, (int)_c.DIO_QUEST_SENS_DC, T.REAL, "dc sensitivity ") },
            { "sens_real", new IP(IF.OPU, (int)_c.DIO_QUEST_SENS_REAL, T.REAL, "dc sens. and real part of ac sensitivity") },
            { "sens_imag", new IP(IF.OPU, (int)_c.DIO_QUEST_SENS_IMAG, T.REAL, "imag part of ac sensitivity ") },
            { "sens_mag", new IP(IF.OPU, (int)_c.DIO_QUEST_SENS_MAG, T.REAL, "sensitivity of ac magnitude") },
            { "sens_ph", new IP(IF.OPU, (int)_c.DIO_QUEST_SENS_PH, T.REAL, "sensitivity of ac phase") },
            { "sens_cplx", new IP(IF.OPU, (int)_c.DIO_QUEST_SENS_CPLX, T.COMPLEX, "ac sensitivity") }
        };

        public bool DIOoff { get; set; } = false;
        public Parameter<double> DIOtemp { get; } = new Parameter<double>(300.15);
        public Parameter<double> DIOinitCond { get; } = new Parameter<double>();
        public Parameter<double> DIOarea { get; } = new Parameter<double>(1.0);
        public int DIOsenParmNo { get; set; } = 0;
        #endregion

        #region Constants
        private const int DIOvoltage = 0;
        private const int DIOcurrent = 1;
        private const int DIOconduct = 2;
        private const int DIOcapCharge = 3;
        private const int DIOcapCurrent = 4;
        private const int DIOsensxp = 5;

        private const int DIORSNOIZ = 0;
        private const int DIOIDNOIZ = 1;
        private const int DIOFLNOIZ = 2;
        private const int DIOTOTNOIZ = 3;
        #endregion

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="name"></param>
        public Diode(string name) : base(name, 2, typeof(DiodeModel)) { }

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="name">The name of the diode</param>
        /// <param name="posnode">The cathode of the diode</param>
        /// <param name="negnode">The anode of the diode</param>
        public Diode(string name, string posnode, string negnode) : this(name)
        {
            Connect(posnode, negnode);
        }

        /// <summary>
        /// Set parameter
        /// </summary>
        protected override void Param(int id, object value, Circuit ckt = null)
        {
            switch ((_c)id)
            {
                case _c.DIO_AREA:
                    DIOarea.Par((double)value);
                    break;
                case _c.DIO_TEMP:
                    DIOtemp.Par((double)value + Circuit.CONSTCtoK);
                    break;
                case _c.DIO_OFF:
                    DIOoff = (bool)value;
                    break;
                case _c.DIO_IC:
                    DIOinitCond.Par((double)value);
                    break;
                case _c.DIO_AREA_SENS:
                    DIOsenParmNo = (int)value;
                    break;
                default:
                    throw new BadParameterException();
            }
        }

        /// <summary>
        /// Ask parameter
        /// </summary>
        protected override object Ask(int id, Circuit ckt = null)
        {
            switch ((_c)id)
            {
                case _c.DIO_OFF:
                    return DIOoff;
                case _c.DIO_IC:
                    return DIOinitCond.Value;
                case _c.DIO_AREA:
                    return DIOarea.Value;
                case _c.DIO_TEMP:
                    return DIOtemp - Circuit.CONSTCtoK;
                case _c.DIO_VOLTAGE:
                    return ckt.State.States[0][DIOstate + DIOvoltage];
                case _c.DIO_CURRENT:
                    return ckt.State.States[0][DIOstate + DIOcurrent];
                case _c.DIO_CAP:
                    return DIOcap;
                case _c.DIO_CHARGE:
                    return ckt.State.States[0][DIOstate + DIOcapCharge];
                case _c.DIO_CAPCUR:
                    return ckt.State.States[0][DIOstate + DIOcapCurrent];
                case _c.DIO_CONDUCT:
                    return ckt.State.States[0][DIOstate + DIOconduct];
                case _c.DIO_POWER:
                    if (ckt.CurrentAnalysis.HasFlag(AnalysisType.DoingAc))
                    {
                        throw new AcCurrentPowerException();
                    }
                    else
                    {
                        return ckt.State.States[0][DIOstate + DIOcurrent] *
                                ckt.State.States[0][DIOstate + DIOvoltage];
                    }
                default:
                    throw new BadParameterException();
            }
        }

        /// <summary>
        /// Set up the diode
        /// </summary>
        /// <param name="ckt"></param>
        public override void Setup(Circuit ckt, Model mod, object p = null)
        {
            DiodeModel model = mod as DiodeModel;

            // Bind nodes
            CircuitNode[] nodes = null;
            if (model.DIOresist == 0.0)
            {
                nodes = BindNodes(ckt);
                DIOposPrimeNode = nodes[0].Number;
            }
            else
            {
                nodes = BindNodes(ckt, NodeType.Voltage);
                DIOposPrimeNode = nodes[2].Number;
            }
            DIOposNode = nodes[0].Number;
            DIOnegNode = nodes[1].Number;

            // Allocate states
            DIOstate = CreateState(ckt, 5);

            // Sensitivity not supported
        }

        /// <summary>
        /// Do temperature dependent calculations
        /// </summary>
        /// <param name="ckt"></param>
        public override void Temperature(Circuit ckt, Model mod, object p = null)
        {
            DiodeModel model = mod as DiodeModel;

            double vte;
            double cbv;
            double xbv;
            double xcbv = 0.0;
            double tol;
            double vt;
            double egfet1, arg1, fact1, pbfact1, pbo, gmaold;
            double fact2, pbfact, arg, egfet, gmanew;

            int iter;
            if (!(p is DiodeModel.DIOTempParameters))
                throw new CircuitException("Invalid parameters specified");
            DiodeModel.DIOTempParameters tp = p as DiodeModel.DIOTempParameters;

            /* loop through all the instances */
            if (!DIOtemp.Given)
                DIOtemp.Value = ckt.Config.Temperature;
            vt = Circuit.CONSTKoverQ * DIOtemp;
            /* this part gets really ugly - I won't even try to
             * explain these equations */
            fact2 = DIOtemp / Circuit.CONSTRefTemp;
            egfet = 1.16 - (7.02e-4 * DIOtemp * DIOtemp) /
                    (DIOtemp + 1108);
            arg = -egfet / (2 * Circuit.CONSTBoltz * DIOtemp) +
                    1.1150877 / (Circuit.CONSTBoltz * (Circuit.CONSTRefTemp + Circuit.CONSTRefTemp));
            pbfact = -2 * vt * (1.5 * Math.Log(fact2) + Circuit.CHARGE * arg);
            egfet1 = 1.16 - (7.02e-4 * model.DIOnomTemp * model.DIOnomTemp) /
                    (model.DIOnomTemp + 1108);
            arg1 = -egfet1 / (Circuit.CONSTBoltz * 2 * model.DIOnomTemp) +
                    1.1150877 / (2 * Circuit.CONSTBoltz * Circuit.CONSTRefTemp);
            fact1 = model.DIOnomTemp / Circuit.CONSTRefTemp;
            pbfact1 = -2 * tp.vtnom * (1.5 * Math.Log(fact1) + Circuit.CHARGE * arg1);
            pbo = (model.DIOjunctionPot - pbfact1) / fact1;
            gmaold = (model.DIOjunctionPot - pbo) / pbo;
            DIOtJctCap = model.DIOjunctionCap /
                    (1 + model.DIOgradingCoeff *
                    (400e-6 * (model.DIOnomTemp - Circuit.CONSTRefTemp) - gmaold));
            DIOtJctPot = pbfact + fact2 * pbo;
            gmanew = (DIOtJctPot - pbo) / pbo;
            DIOtJctCap *= 1 + model.DIOgradingCoeff *
                    (400e-6 * (DIOtemp - Circuit.CONSTRefTemp) - gmanew);

            DIOtSatCur = model.DIOsatCur * Math.Exp(
                    ((DIOtemp / model.DIOnomTemp) - 1) *
                    model.DIOactivationEnergy / (model.DIOemissionCoeff * vt) +
                    model.DIOsaturationCurrentExp / model.DIOemissionCoeff *
                    Math.Log(DIOtemp / model.DIOnomTemp));
            /* the defintion of f1, just recompute after temperature adjusting
             * all the variables used in it */
            DIOtF1 = DIOtJctPot *
                    (1 - Math.Exp((1 - model.DIOgradingCoeff) * tp.xfc)) /
                    (1 - model.DIOgradingCoeff);
            /* same for Depletion Capacitance */
            DIOtDepCap = model.DIOdepletionCapCoeff *
                    DIOtJctPot;
            /* and Vcrit */
            vte = model.DIOemissionCoeff * vt;
            DIOtVcrit = vte * Math.Log(vte / (Circuit.CONSTroot2 * DIOtSatCur));
            /* and now to copute the breakdown voltage, again, using
             * temperature adjusted basic parameters */
            if (model.DIObreakdownVoltage.Given)
            {
                cbv = model.DIObreakdownCurrent;
                if (cbv < DIOtSatCur * model.DIObreakdownVoltage / vt)
                {
                    cbv = DIOtSatCur * model.DIObreakdownVoltage / vt;
                    CircuitWarning.Warning(this, string.Format("Diode {0}, model {1}: Breakdown current increased to {2} to resolve incompatibility with specified saturation current", Name, model.Name, cbv));
                    xbv = model.DIObreakdownVoltage;
                }
                else
                {
                    tol = ckt.Config.RelTol * cbv;
                    xbv = model.DIObreakdownVoltage - vt * Math.Log(1 + cbv /
                            DIOtSatCur);
                    iter = 0;
                    for (iter = 0; iter < 25; iter++)
                    {
                        xbv = model.DIObreakdownVoltage - vt * Math.Log(cbv /
                                DIOtSatCur + 1 - xbv / vt);
                        xcbv = DIOtSatCur * (Math.Exp((model.DIObreakdownVoltage
                                - xbv) / vt) - 1 + xbv / vt);
                        if (Math.Abs(xcbv - cbv) <= tol) goto matched;
                    }
                    CircuitWarning.Warning(this, string.Format("Diode {0}, model {1}: Unable to match forward and reverse diode regions: bv = {2}, ibv = {3}", Name, model.Name, xbv, xcbv));
                }
            matched:
                DIOtBrkdwnV = xbv;
            }
        }

        /// <summary>
        /// Load the circuit
        /// </summary>
        /// <param name="ckt"></param>
        public override void Load(Circuit ckt, Model mod, object p = null)
        {
            DiodeModel model = mod as DiodeModel;

            double arg;
            double capd;
            double cd;
            double cdeq;
            double cdhat;
            // double ceq;
            double csat;    /* area-scaled saturation current */
            double czero;
            double czof2;
            double delvd;   /* change in diode voltage temporary */
            double evd;
            double evrev;
            double gd;
            // double geq;
            double gspr;    /* area-scaled conductance */
            double sarg;
            // double tol;     /* temporary for tolerence calculations */
            double vd;      /* current diode voltage */
            double vdtemp;
            double vt;      /* K t / Q */
            double vte;
            bool Check;
            // bool SendCond = false;

            CircuitState state = ckt.State;

            /*
             *     this routine loads diodes for dc and transient analyses.
             */


            /* if(ckt->CKTsenInfo){
                if((ckt->CKTsenInfo->SENstatus == PERTURBATION)
                        && (DIOsenPertFlag == OFF))continue;
                SenCond = DIOsenPertFlag;
            } */

            csat = DIOtSatCur * DIOarea;
            gspr = model.DIOconductance * DIOarea;
            vt = Circuit.CONSTKoverQ * DIOtemp;
            vte = model.DIOemissionCoeff * vt;
            /*  
             *   initialization 
             */

            /* if(SenCond){
                if((ckt->CKTsenInfo->SENmode == TRANSEN)&&
                        ((ckt.Mode & Circuit.Modes.InitTran) != 0)) {
                    vd = state.States[1][DIOstate + DIOvoltage];
                } else{
                    vd = state.States[0][DIOstate + DIOvoltage];
                }
                goto next1;
            } */

            Check = true;
            if ((ckt.Mode & Circuit.Modes.InitSmSig) != 0)
            {
                vd = state.States[0][DIOstate + DIOvoltage];
            }
            else if ((ckt.Mode & Circuit.Modes.InitTran) != 0)
            {
                vd = state.States[1][DIOstate + DIOvoltage];
            }
            else if (((ckt.Mode & Circuit.Modes.InitJct) != 0) &&
                  ((ckt.Mode & Circuit.Modes.TranOp) != 0) && ((ckt.Mode & Circuit.Modes.Uic) != 0))
            {
                vd = DIOinitCond;
            }
            else if (((ckt.Mode & Circuit.Modes.InitJct) != 0) && DIOoff)
            {
                vd = 0;
            }
            else if ((ckt.Mode & Circuit.Modes.InitJct) != 0)
            {
                vd = DIOtVcrit;
            }
            else if (((ckt.Mode & Circuit.Modes.InitFix) != 0) && DIOoff)
            {
                vd = 0;
            }
            else
            {
                vd = state.Solution[DIOposPrimeNode] -
                        state.Solution[DIOnegNode];
                delvd = vd - state.States[0][DIOstate + DIOvoltage];
                cdhat = state.States[0][DIOstate + DIOcurrent] +
                        state.States[0][DIOstate + DIOconduct] * delvd;
                /*  
                 *   bypass if solution has not changed
                 */

                /*
                 *   limit new junction voltage
                 */
                if ((model.DIObreakdownVoltage.Given) &&
                        (vd < Math.Min(0, -DIOtBrkdwnV + 10 * vte)))
                {
                    vdtemp = -(vd + DIOtBrkdwnV);
                    vdtemp = pnjlim(vdtemp,
                            -(state.States[0][DIOstate + DIOvoltage] +
                            DIOtBrkdwnV), vte,
                            DIOtVcrit, ref Check);
                    vd = -(vdtemp + DIOtBrkdwnV);
                }
                else
                {
                    vd = pnjlim(vd, state.States[0][DIOstate + DIOvoltage],
                            vte, DIOtVcrit, ref Check);
                }
            }
            /*
             *   compute dc current and derivatives
             */
            // next1: 
            if (vd >= -3 * vte)
            {
                evd = Math.Exp(vd / vte);
                cd = csat * (evd - 1) + ckt.Config.Gmin * vd;
                gd = csat * evd / vte + ckt.Config.Gmin;
            }
            else if ((DIOtBrkdwnV == 0.0) || vd >= -DIOtBrkdwnV)
            {
                arg = 3 * vte / (vd * Circuit.CONSTE);
                arg = arg * arg * arg;
                cd = -csat * (1 + arg) + ckt.Config.Gmin * vd;
                gd = csat * 3 * arg / vd + ckt.Config.Gmin;
            }
            else
            {
                evrev = Math.Exp(-(DIOtBrkdwnV + vd) / vte);
                cd = -csat * evrev + ckt.Config.Gmin * vd;
                gd = csat * evrev / vte + ckt.Config.Gmin;
            }
            if (((ckt.Mode & (Circuit.Modes.Tran | Circuit.Modes.Ac | Circuit.Modes.InitSmSig)) != 0) ||
                    ((ckt.Mode & Circuit.Modes.TranOp) != 0) && ((ckt.Mode & Circuit.Modes.Uic) != 0))
            {
                /*
                 *   charge storage elements
                 */
                czero = DIOtJctCap * DIOarea;
                if (vd < DIOtDepCap)
                {
                    arg = 1 - vd / model.DIOjunctionPot;
                    sarg = Math.Exp(-model.DIOgradingCoeff * Math.Log(arg));
                    state.States[0][DIOstate + DIOcapCharge] =
                            model.DIOtransitTime * cd + model.DIOjunctionPot *
                            czero * (1 - arg * sarg) / (1 - model.DIOgradingCoeff);
                    capd = model.DIOtransitTime * gd + czero * sarg;
                }
                else
                {
                    czof2 = czero / model.DIOf2;
                    state.States[0][DIOstate + DIOcapCharge] =
                            model.DIOtransitTime * cd + czero * DIOtF1 + czof2 *
                            (model.DIOf3 * (vd - DIOtDepCap) +
                            (model.DIOgradingCoeff / (model.DIOjunctionPot +
                            model.DIOjunctionPot)) * (vd * vd - DIOtDepCap *
                            DIOtDepCap));
                    capd = model.DIOtransitTime * gd + czof2 * (model.DIOf3 +
                            model.DIOgradingCoeff * vd / model.DIOjunctionPot);
                }
                DIOcap = capd;

                /*
                 *   store small-signal parameters
                 */
                if ((!((ckt.Mode & Circuit.Modes.TranOp) != 0)) ||
                        (!((ckt.Mode & Circuit.Modes.Uic) != 0)))
                {
                    if ((ckt.Mode & Circuit.Modes.InitSmSig) != 0)
                    {
                        state.States[0][DIOstate + DIOcapCurrent] = capd;

                        /* if(SenCond){
                            state.States[0][DIOstate + DIOcurrent] = cd;
                            state.States[0][DIOstate + DIOconduct] = gd;
                        } */
                        return;
                    }

                    /*
                     *   transient analysis
                     */
                    /* if(SenCond && (ckt->CKTsenInfo->SENmode == TRANSEN)){
                        state.States[0][DIOstate + DIOcurrent] = cd;
                        return;
                    } */

                    if ((ckt.Mode & Circuit.Modes.InitTran) != 0)
                    {
                        state.States[1][DIOstate + DIOcapCharge] =
                                state.States[0][DIOstate + DIOcapCharge];
                    }
                    var result = ckt.Integration.Integrate(state, DIOstate + DIOcapCharge, capd);
                    gd = gd + result.Geq;
                    cd = cd + state.States[0][DIOstate + DIOcapCurrent];
                    if ((ckt.Mode & Circuit.Modes.InitTran) != 0)
                    {
                        state.States[1][DIOstate + DIOcapCurrent] =
                                state.States[0][DIOstate + DIOcapCurrent];
                    }
                }
            }

            // if(SenCond) goto next2;

            /*
             *   check convergence
             */
            if ((!((ckt.Mode & Circuit.Modes.InitFix) != 0)) || (!(DIOoff)))
            {
                if (Check)
                    ckt.IsCon = false;
            }
            // next2:      
            state.States[0][DIOstate + DIOvoltage] = vd;
            state.States[0][DIOstate + DIOcurrent] = cd;
            state.States[0][DIOstate + DIOconduct] = gd;

            // if(SenCond)  continue;


            // load:

            /*
             *   load current vector
             */
            cdeq = cd - gd * vd;
            state.Rhs[DIOnegNode] += cdeq;
            state.Rhs[DIOposPrimeNode] -= cdeq;

            /*
             *   load matrix
             */
            state.Matrix[DIOposNode, DIOposNode] += gspr;
            state.Matrix[DIOnegNode, DIOnegNode] += gd;
            state.Matrix[DIOposPrimeNode, DIOposPrimeNode] += (gd + gspr);
            state.Matrix[DIOposNode, DIOposPrimeNode] -= gspr;
            state.Matrix[DIOnegNode, DIOposPrimeNode] -= gd;
            state.Matrix[DIOposPrimeNode, DIOposNode] -= gspr;
            state.Matrix[DIOposPrimeNode, DIOnegNode] -= gd;
        }

        /// <summary>
        /// Load the diode for AC simulation
        /// </summary>
        /// <param name="ckt"></param>
        public override void AcLoad(Circuit ckt, Model mod, object p = null)
        {
            DiodeModel model = mod as DiodeModel;

            double gspr;
            double geq;
            double xceq;
            CircuitStateComplex state = ckt.StateComplex;
            CircuitState cstate = ckt.State;

            gspr = model.DIOconductance * DIOarea; // [EDIT] Changed model.DIOresist to model.DIOconductance
            geq = cstate.States[0][DIOstate + DIOconduct];
            xceq = cstate.States[0][DIOstate + DIOcapCurrent];

            state.Matrix[DIOposNode, DIOposNode] += gspr;
            // *(DIOposPosPtr) += gspr;
            state.Matrix[DIOnegNode, DIOnegNode] += geq + xceq * state.Omega;
            // *(DIOnegNegPtr) += geq + xceq * s->real;
            // *(DIOnegNegPtr + 1) += xceq * s->imag;
            state.Matrix[DIOposPrimeNode, DIOposPrimeNode] += geq + gspr + xceq * state.Omega;
            // *(DIOposPrimePosPrimePtr) += geq + gspr + xceq * s->real;
            // *(DIOposPrimePosPrimePtr + 1) += xceq * s->imag;
            state.Matrix[DIOposNode, DIOposPrimeNode] -= gspr;
            // *(DIOposPosPrimePtr) -= gspr;
            state.Matrix[DIOnegNode, DIOposPrimeNode] -= geq + xceq * state.Omega;
            // *(DIOnegPosPrimePtr) -= geq + xceq * s->real;
            // *(DIOnegPosPrimePtr + 1) -= xceq * s->imag;
            state.Matrix[DIOposPrimeNode, DIOposNode] -= gspr;
            // *(DIOposPrimePosPtr) -= gspr;
            state.Matrix[DIOposPrimeNode, DIOnegNode] -= geq + xceq * state.Omega;
            // *(DIOposPrimeNegPtr) -= geq + xceq * s->real;
            // *(DIOposPrimeNegPtr + 1) -= xceq * s->imag;
        }

        /// <summary>
        /// Load noise contributions
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="mode"></param>
        /// <param name="operation"></param>
        /// <param name="data"></param>
        /// <param name="OnDens"></param>
        /// <param name="p"></param>
        public override void NoiseLoad(Circuit ckt,
            Noise.OutputMode mode,
            Noise.Operation operation,
            NoiseData data,
            ref double OnDens,
            Model mod, object p = null)
        {
            DiodeModel model = mod as DiodeModel;

            double tempOnoise;
            double tempInoise;
            double[] noizDens = new double[DIOnVar.Length];
            double[] lnNdens = new double[DIOnVar.Length];
            double tmp = 0.0;
            // int error;
            int i;
            NoiseJob curjob = ckt.CurrentJob as NoiseJob;
            CircuitState cstate = ckt.State;
            CircuitStateComplex state = ckt.StateComplex;

            /* define the names of the noise sources */

            string[] DIOnNames = new string[] {       /* Note that we have to keep the order */
	            "_rs",              /* noise due to rs */        /* consistent with the index definitions */
	            "_id",              /* noise due to id */        /* in DIOdefs.h */
	            "_1overf",          /* flicker (1/f) noise */
	            ""                  /* total diode noise */
            };

            switch (operation)
            {

                case Noise.Operation.N_OPEN:

                    /* see if we have to to produce a summary report */
                    /* if so, name all the noise generators */

                    if (curjob.NStpsSm != 0)
                    {
                        switch (mode)
                        {

                            case Noise.OutputMode.N_DENS:
                                for (i = 0; i < DIOnVar.Length; i++)
                                    data.namelist.Add(String.Format("onoise_{0}{1}", Name, DIOnNames[i]));
                                break;

                            case Noise.OutputMode.INT_NOIZ:
                                for (i = 0; i < DIOnVar.Length; i++)
                                {
                                    data.namelist.Add(String.Format("onoise_total_{0}{1}", Name, DIOnNames[i]));
                                    data.namelist.Add(String.Format("inoise_total_{0}{1}", Name, DIOnNames[i]));
                                }
                                break;
                        }
                    }
                    break;

                case Noise.Operation.N_CALC:
                    switch (mode)
                    {

                        case Noise.OutputMode.N_DENS:
                            NevalSrc(ckt, ref noizDens[DIORSNOIZ], ref lnNdens[DIORSNOIZ],
                                 Noise.Sources.THERMNOISE, DIOposPrimeNode, DIOposNode,
                                 model.DIOconductance * DIOarea);

                            NevalSrc(ckt, ref noizDens[DIOIDNOIZ], ref lnNdens[DIOIDNOIZ],
                                     Noise.Sources.SHOTNOISE, DIOposPrimeNode, DIOnegNode,
                                 cstate.States[0][DIOstate + DIOcurrent]);

                            NevalSrc(ckt, ref noizDens[DIOFLNOIZ], ref tmp,
                                 Noise.Sources.N_GAIN, DIOposPrimeNode, DIOnegNode,
                                 (double)0.0);
                            noizDens[DIOFLNOIZ] *= model.DIOfNcoef *
                                 Math.Exp(model.DIOfNexp *
                                 Math.Log(Math.Max(Math.Abs(cstate.States[0][DIOstate + DIOcurrent]), Noise.N_MINLOG))) /
                                 data.freq;
                            lnNdens[DIOFLNOIZ] =
                                 Math.Log(Math.Max(noizDens[DIOFLNOIZ], Noise.N_MINLOG));

                            noizDens[DIOTOTNOIZ] = noizDens[DIORSNOIZ] +
                                            noizDens[DIOIDNOIZ] +
                                            noizDens[DIOFLNOIZ];
                            lnNdens[DIOTOTNOIZ] =
                                 Math.Log(Math.Max(noizDens[DIOTOTNOIZ], Noise.N_MINLOG));

                            OnDens += noizDens[DIOTOTNOIZ];

                            if (data.delFreq == 0.0)
                            {

                                /* if we haven't done any previous integration, we need to */
                                /* initialize our "history" variables                      */

                                for (i = 0; i < DIOnVar.Length; i++)
                                {
                                    DIOnVar[i].LNLSTDENS = lnNdens[i];
                                }

                                /* clear out our integration variables if it's the first pass */

                                if (data.freq == curjob.NstartFreq)
                                {
                                    for (i = 0; i < DIOnVar.Length; i++)
                                    {
                                        DIOnVar[i].OUTNOIZ = 0.0;
                                        DIOnVar[i].INNOIZ = 0.0;
                                    }
                                }
                            }
                            else
                            {   /* data.delFreq != 0.0 (we have to integrate) */

                                /* To insure accurracy, we have to integrate each component separately */

                                for (i = 0; i < DIOnVar.Length; i++)
                                {
                                    if (i != DIOTOTNOIZ)
                                    {
                                        tempOnoise = NoiseIntegration.Integrate(noizDens[i], lnNdens[i],
                                              DIOnVar[i].LNLSTDENS, data);
                                        tempInoise = NoiseIntegration.Integrate(noizDens[i] * data.GainSqInv,
                                              lnNdens[i] + data.lnGainInv,
                                              DIOnVar[i].LNLSTDENS + data.lnGainInv,
                                              data);
                                        DIOnVar[i].LNLSTDENS = lnNdens[i];
                                        data.outNoiz += tempOnoise;
                                        data.inNoise += tempInoise;
                                        if (curjob.NStpsSm != 0)
                                        {
                                            DIOnVar[i].OUTNOIZ += tempOnoise;
                                            DIOnVar[DIOTOTNOIZ].OUTNOIZ += tempOnoise;
                                            DIOnVar[i].INNOIZ += tempInoise;
                                            DIOnVar[DIOTOTNOIZ].INNOIZ += tempInoise;
                                        }
                                    }
                                }
                            }
                            if (data.prtSummary)
                            {
                                for (i = 0; i < DIOnVar.Length; i++)
                                {     /* print a summary report */
                                    data.outpVector[data.outNumber++] = noizDens[i];
                                }
                            }
                            break;

                        case Noise.OutputMode.INT_NOIZ:        /* already calculated, just output */
                            if (curjob.NStpsSm != 0)
                            {
                                for (i = 0; i < DIOnVar.Length; i++)
                                {
                                    data.outpVector[data.outNumber++] = DIOnVar[i].OUTNOIZ;
                                    data.outpVector[data.outNumber++] = DIOnVar[i].INNOIZ;
                                }
                            }    /* if */
                            break;
                    }    /* switch (mode) */
                    break;

                case Noise.Operation.N_CLOSE:
                    return;         /* do nothing, the main calling routine will close */
            }    /* switch (operation) */
        }
    }
}
